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I. INTRODUCTION 

The very first suggestion to study J / ip production and 
suppression in heavy-ion collisions [1] was based on a two- 
fold observation: First, the energy density accumulated 
in such collisions is expected to be high enough that the 
medium formed in the collision presumably undergoes a 
transition from hadronic matter to a quark-gluon plasma 
(QGP). Second, J/ijj is bound essentially by the con- 
finement force which, when screened in the QGP, loses 
its binding character. Produced J/ips dissolve and their 
charm quarks end up in pairs of D mesons. The expecta- 
tion is therefore that the J / ip production cross section is 
suppressed in a hot and dense environment as compared 
with extrapolations from proton-proton (jyp) collisions. 

Indeed, lattice simulations [2] indicate that QCD with 
(2+1) flavors undergoes a transition at a temperature of 
about 170 MeV [3] from a confined hadronic phase to a 
phase where quarks and gluons constitute the relevant 
active degrees of freedom. The possible existence of this 
new state of matter has raised enormous interest, and the 
search for experimental evidence of its creation continues 
since about 20 years. Yet there is still no unambiguous 
proof that the QGP has been produced in the laboratory. 
A primary reason for this is that any QGP signature, in- 
cluding J/ip suppression, is folded with the time evolu- 
tion of the fireball created in such a collision. Moreover, 
the evolution continues after the supposedly produced 
plasma undergoes a transition from partonic to hadronic 
degrees of freedom. Any information on the early stage is 
then hidden behind signals from the hadronic phase. Fi- 
nally, the most important uncertainty is perhaps the poor 
knowledge of the process of thermalization. The partic- 
ularly interesting case of J/ip suppression is thus subject 
of a vigorous experimental search and theoretical debate. 
The first systematic signs of suppression [4] were later on 
explained by more conventional mechanisms [5], already 
present in proton-nucleus (pA) reactions. Only with the 
beginning of the experimental search using Pb+Pb colli- 
sions the first signals of anomalous suppression, beyond 
extrapolations from pA collisions, were seen [6]. 

While much phenomenology has been developed in or- 



der to explain the latest SPS data [7], it is still debated 
whether or not anomalous suppression can be attributed 
to QGP formation. In fact, realistic models can be con- 
structed which incorporate well known nuclear effects 
such as initial state gluon radiation [8], color excitation 
[9], initial state parton energy loss [10] and coherence ef- 
fects [11]. All of these can potentially account for the 
observed data because they provide the necessary non- 
linear dependence on the number of participating parti- 
cles which becomes significant only when going beyond 
pA collisions. 

In view of the uncertainties inherent in the construc- 
tion of models which consider such effects, we attempt an 
approach to the problem which maximizes the number of 
independent physical constraints. The usual strategy is 
to attribute the missing suppression to medium effects, 
exclusively designed to handle the J/ip case alone. For 
this procedure to be conclusive one needs a very reli- 
able baseline for the description of the production pro- 
cess. Here we adopt a different and novel perspective and 
make use of the knowledge of the medium evolution as 
inferred by a variety of other observables. We then ex- 
plore whether the same description is consistent with the 
observed J/ip measurements. We stress that this is pos- 
sible because the specific time evolution of the medium is 
not constructed or fitted in order to reproduce the J/ip 
suppression effect, but it is constrained beforehand by 
other independent data. 

We treat charmonium production in nucleus-nucleus 
(AB) collisions as a two-step process, factorizing direct 
production from the subsequent evolution in the medium. 
The first part includes the conventional description [13] 
of nuclear effects within the Glauber model framework, 
with values of the absorption cross section either extrap- 
olated from the suppression observed experimentally al- 
ready in pA collisions, or fitted to reproduce the results 
of more sophisticated computations [11, 12]. The sec- 
ond part is a description within kinetic theory of the 
interaction of J / ips with the different degrees of freedom 
that populate the evolving medium at different times. In 
this framework interactions are realized by cross sections, 
while the evolution of the medium itself is constrained 
by an independent description of hadronic spectra. With 
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this procedure we are able to specify the nature of the 
medium and to follow the detailed time evolution of J / ip 
mesons. 

It might seem natural at this stage to also incorporate 
the mentioned Debye screening of the inter-quark poten- 
tial, for example by assuming complete J/ip dissolution 
above a certain temperature [14, 15]. However, it is not 
clear what an all-embracing description of the interaction 
of the J/ip with the medium should be. Screening arises 
from the interaction of the virtual gluons that bind the 
J/ip with the gluons of the environment. Such a process 
can also be viewed as scattering of the J/ip, that fluc- 
tuates into a cc pair and non-perturbative multi-gluon 
exchanges, off the thermalized gluons [16]. Taking into 
account both screening and kinetic scattering can amount 
to a certain degree to double counting. Disentangling 
these effects even qualitatively, however, is a highly non- 
trivial task which requires a dedicated analysis. At this 
stage, we can say that screening, as inferred from lat- 
tice QCD, is a purely static concept. Since J/ips are not 
produced at rest with respect to the medium, a kinetic 
approach seems to be appropriate within the present con- 
text. Obviously, a more rigorous treatment of this ques- 
tion is required in the future. 

The paper is structured in four sections. First, the ini- 
tial conditions for the solution of the kinetic equation are 
discussed. We provide a simple extension of the parton 
model prescription for the pp — > ccX process, scaling 
the open charm production cross section with the aver- 
age number of collisions. At the same time charmonium 
production is computed treating nuclear suppression as 
effective Glauber absorption. We take into account the 
role of excited charmonium states ip' and \ which con- 
tribute sequentially to J/ip production. We will often 
denote charmonia generically by as a shorthand nota- 
tion. 

In the following section a thorough description of the 
produced medium is given. This medium is assumed to 
thermalize and then described within a thermodynami- 
cally consistent quasi-particle picture [17], incorporating 
important aspects of confinement. This approach treats 
quarks and gluons as massive thermal quasi-particles, 
with their properties determined to reproduce lattice 
QCD results. The driving force of the transition, the 
confinement/dcconfincment process, is given a statisti- 
cal meaning in terms of a reduction of the thermally ac- 
tive partonic degrees of freedom as the critical tempera- 
ture is approached from above. The fireball evolution is 
fixed by requiring agreement with hadronic observables 
at freeze-out [18]. Thermodynamics is calculated along 
volumes of constant proper time under the assumption 
of total entropy conservation and using the equation of 
state (EoS) from lattice QCD as described by the quasi- 
particle model. In a self-consistent calculational proce- 
dure, we determine temperature and volume of the sys- 
tem as a function of proper time. 

The main part of the paper is then focused on the ki- 
netic description of charmonium evolution in the fireball. 



A transport equation is presented, including dissociation 
processes in the collision term. In the QGP, <J/ dissocia- 
tion is provided by the well known Bhanot-Peskin cross 
section [19, 20]. We also take into account the possibil- 
ity of ^ regeneration in the QGP. As recently suggested 
by several authors [21-24], this process can potentially 
overwhelm dissociation, leading to enhancement once 
the heavy ion collision energy is large enough. The coa- 
lescence process of cc quark pairs into <!/ is controlled by 
the cross section obtained from that of dissociation us- 
ing detailed balance. A simplified, approximate kinetic 
equation is then obtained. Its solution can be written in 
closed form, with time dependence entering through the 
evolution of fireball temperature and volume. 

Finally, results for * production at SPS and RHIC 
are presented and a thorough discussion of theoretical 
uncertainties concludes the paper. 

II. CHARM PRODUCTION OFF NUCLEI: 
INITIAL CONDITIONS 

A. Open charm 

The production of charmed quarks is commonly de- 
scribed within pcrturbative QCD. The pcrturbative ap- 
proach is strictly valid only for processes involving large 
gluon virtualities, usually provided by high-momentum 
exchange or large masses. In the present case the small- 
est scale is that of the charm quark mass m c ~ 1.5 GeV, 
implying that perturbation theory at lowest order re- 
ceives important corrections. Moreover, computations 
are accurate only at large enough transverse momenta 
Pt, whereas small values of pt contribute significantly to 
total yields. 

In the following discussion we restrict ourselves to the 
leading order treatment presented in [25, 26], with suit- 
able adjustments in order to meet phenomenology. We 
consider the leading processes qq — ► cc and gg — > cc. 
In terms of their elementary differential cross sections 
d&i/dt, the spectrum of c quarks produced in pp colli- 
sions at rapidity y c and transverse momentum is 

1 i—u,d,s,g 

where x\.2 — (m^/y/s) [cxp(±y e ) + cxp(±7/ c )] are the mo- 
mentum fractions of the partons in the colliding protons 
and Fi/ p = Xifi/ p . The rapidity variable is defined as 
usual as y = 1/2 log[ (E + p z )/(E — p z ) ]. The factoriza- 
tion scale is taken at /i c = 1.4 GeV, which is of the order 
of the c quark mass, and a scaling factor K — 2 is used. 
We employ the GRV94LO [27] parton distributions f^ p 
and neglect the effect of intrinsic transverse parton mo- 
mentum. Moreover we set in cq. (1), where appropriate, 
m c = [i c . In this way it is possible to effectively repro- 
duce to good accuracy the next-to-leading order compu- 
tations for total (y c - and pT-integrated) charm produc- 
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tion as given in [26] and conveniently parametrized in 
[24] as 



<ifc(*) = a° c (l m /Vi) 8 - 185 (^/m ) 



1.132 



(2) 



where cr° = 3.392 and m = 2.984 GeV. 

To estimate the spectrum of c quarks in AB collisions, 
the simplest approach is to scale the pp result with the 
overlap function 



Tab (b) = j d 2 r T A {r)T B (r) 



(3) 



where Ta,b{V) = J dz pa,b{z,V) are the usual overlap 
functions, expressed in terms of longitudinal integrals 
over the nuclear densities Pa,b, and r — \b — r\, being 
f the transverse coordinate and b the impact parameter. 

On the other hand, as the collision energy increases, 
shadowing effects are expected to become important, re- 
ducing the total yield [28]. Another correction to be in- 
troduced is the broadening effect on the colliding partons' 
transverse momenta, resulting in a broader pt spectrum 
of produced charmed quarks. For simplicity we neglect 
these effects, thereby obtaining an upper limit for the 
c quark production cross section. The charmed quark 
spectrum in AB collisions is then computed as 



dN AB /ja = d(J PP 



dy c dp\ 



dy c dp\ 



T AB (b). 



(4) 



The total number of charm quarks produced, for exam- 
ple, in a Au+Au collision can be obtained by integrating 
the spectrum given in the latter equation. The result is 
shown in Fig. (1). One observes a great difference be- 
tween the values at SPS and RH1C energies, amounting 
to nearly two orders of magnitude. Moreover, note that 

N c AB (s) = a pp ^ c (s) T AB (b) = NUs) N%° B \b) , (5) 



where Ng p 



^ pp 



AB (b 



«PP^c/v p n p and N%° B l (b) = a™ T AB , being 
the inelastic pp cross section. Since for large nuclei 
0) > 1000, it is clear how many more c quarks 
are produced in a AB collision with respect to the pp 
case. 



B. Hidden Charm 

The description of charmonium production in nuclear 
collisions is a more complicated task as compared to 
open cc production. At the root of the problem lies the 
fact that exclusive production of composite particles in 
hadronic collisions is basically a non-perturbative pro- 
cess. Only at large transverse momenta it is possible to 
make robust predictions for the spectra. Nevertheless, 
at least at the phenomenological level, much work has 
been done in order to understand the physics underlying 
the results of several experiments. Let us examine pp 
collisions 
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FIG. 1: Number of cc pairs produced in a Au + Au collision as 
function of the center of mass energy, computed by integrat- 
ing the spectrum given in eq. (4). Curves compare the NLO 
parton model result (dotted curve) with modified LO calcula- 
tions with n c = 1.4 GeV (solid). From SPS to RHIC energy 
the number of produced pairs grows by nearly two orders of 
magnitude 



We consider first of all the total production cross sec- 
tion and make use of the parametrization given in [29]. 
Since it was obtained by fitting low-energy data (y/s < 50 
GeV), we modified it to simulate the high energy behav- 
ior in terms of NLO effects in the color evaporation model 
[30]. We parametrize 



0) = °j/</> (! - mj/4,/Vs) 9 (1 + 0.6 y/s/mj/j,) , 



with 



'J/4, 



100 nb and 



3.1 GeV. 



(6) 

The above 



formula gives essentially identical results as that of [29] 
in the low energy region and is in line with the re- 
cent PHENIX measurement at RHIC [31]. The rapidity 
modulation can be inferred from the relation da/dy ~ 
xig(xi) xig[x-i) where g{x) ~ (1 — x) 5 /x is taken to 
be the gluon distribution in the proton and x\^ — 
(m-q, j \fs) exp(±y). We can then write the transverse- 
momentum-integrated ^ production spectrum in pp col- 
lisions as 



da. 



pp^J/iP 
dy 



= Pp P ^j/i,{s) F(s,y) . 



(7) 



where the y-dependent part is 
F(s,y) = C(s) (1 - Xl ) 5 (1 - x 2 f 

r i 5 

= C{s) 1 - 2(m J/l/ ,/ v / i)coshj/ + m J/lf ,/s^ (8) 

while C(s) is chosen to satisfy the normalization con- 
straint Jdy F(s,y) = 1. The transverse momentum de- 
pendence will not be needed in the following, except for 
the value of (jij) which can be taken directly from ex- 
periment. 
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We now consider nuclear effects, starting with the sim- 
pler case of proton- nucleus (pA) collisions. It has been 
shown that the experimental results on J/ip production 
can be described using 

o v a^ji^ = J d 2 b T A (b) S a A bs (b) (9) 

for the total production cross section. The factor 



nao. 



(b) 



1 — exp 



-a 



abs 

J/i/jN 



T A (b) 



T abs 

' j/<4>n 



T A (b) 



(10) 



is the survival probability for J/ip to escape the nucleus 
without being dissociated. It includes the effective ab- 
sorption cross section (Tj b , s ^ N , a quantity of the order of 3 
mb for mid-rapidity J ftps as measured at Ei a b — 800 GeV 
at Fermilab, while it amounts to 4— 6 mb for mid-rapidity 
J/ips as measured at E lab = 200 - 400 GeV at the SPS. 
The absorption cross section parametrizes various poorly 
known effects, with varying importance depending on the 
collision energy. Among those effects are the presence of 
color non-singlet degrees of freedom in the dynamics of 
colliding nucleons, initial state parton energy loss and co- 
herence length and shadowing effects. Moreover, u°j bs ^ N 
contains information on the nuclear absorption of the ex- 
cited states ip' and \c which feed into the observed J/ip 
yield. 

A common property of all of the above effects is the 
approximate linear dependence on the length of the path 
which J/ip traverses in the nuclei. It is therefore jus- 
tified, to first approximation, to use the same "optical" 
absorption formula eq. (10), provided a re-scaling and re- 
intcrprctation of af^ N -> a»™ and Sff B - S% U B C is 
done. 

When looking at J/ip production in AB collisions, one 
can estimate the cross section for a given impact param- 
eter by generalizing eq. (9). The effects of the produced 
medium are the central topic of this paper and will be 
thoroughly discussed in the following sections. Neglect- 
ing them for the moment, one obtains 



dN 



j/4> 

AB 



dy 



(b) 



da 



dy 



T AB (b)S^E C (b), (11) 



where all nuclear effects are included in the suppression 
function 



S^E C (b)=T AB (b) J d 2 r T A (r) Sr C (r) T B (f) S™ c (r) 

(12) 

which has the obvious property S^ B C < 1 and S^ B C — > 

i if NUC n 

Since nuclear effects depend on energy, we have chosen 
(?j/^ C N {so) = 4.5 mb at the SPS energy = 17.3 GeV 
{Eiab = 158 GeV) in order to be in agreement with the 
pA measurement, and assumed the relation 

^%(s) = ^ C N (s ) (s/ So ) x (13) 



with A = 0.12 in order to simulate nuclear effects as pre- 
dicted in [11]. 



III. THERMODYNAMICS OF THE PRODUCED 
MEDIUM 

A. Quasi-particle description of the QGP and the 
Hadronic EoS 

Before we address the question of how to describe the 
medium produced by a nuclear collision, we need to iden- 
tify its relevant degrees of freedom. On the theoretical 
side, information about the QCD medium comes from 
lattice calculations at finite temperature. Simulations of 
the pressure, entropy and energy density [2] show a dras- 
tic increase in the number of degrees of freedom at a 
critical temperature of about 170 MeV (for two massless 
flavors [3]) to about 80% of the value of an ideal gas of 
quarks and gluons. 

A perturbative description of the QGP is, however fu- 
tile, since thermal perturbation theory is in general in- 
sufficient for all temperatures of interest. This is evi- 
dent, for example, from calculations of the free energy of 
the QGP [32] or the photon self-energy in the thermal 
medium [33]. Furthermore, close to T c we expect intrin- 
sically non-perturbative dynamics to enter: the confinc- 
ment/deconfinement transition and spontaneous chiral 
symmetry breaking are not accountable for in an expan- 
sion in the coupling constant. In view of these facts, we 
use a more phenomenological approach to QCD thermo- 
dynamics which nevertheless goes considerably beyond 
commonly used ideal gas models. We have shown re- 
cently that it is possible to describe the EoS of hot 
QCD to a very good approximation by the EoS of a 
gas of quasi-particles with thermally generated masses, 
incorporating confinement effectively by a temperature- 
dependent, reduced number of active degrees of free- 
dom. This reduction is caused by the formation of heavy 
hadrons or glueballs, with large masses compared to T c , 
and accounts for the fact that the QGP does not resemble 
at all an ideal gas of quarks and gluons as the hadroniza- 
tion temperature is approached. Here we give a short 
summary of the method and refer the reader to [17] for 
a more detailed discussion. 

From asymptotic freedom, we expect that at ex- 
tremely high temperatures the plasma consists of quasi- 
free quarks and gluons. Hard thermal loop (HTL) pertur- 
bative calculations find, for thermal momenta, spectral 
functions of the form S{E 2 -k 2 -m 2 {T)) with m(T) - gT 
[34] . As long as the spectral functions of the thermal exci- 
tations at lower temperatures resemble qualitatively this 
asymptotic form, a quasi-particle description is expected 
to be applicable. QCD dynamics is then incorporated 
in the thermal masses of the effective quarks and glu- 
ons, plus an extra function, required by thermodynamical 
consistency, which plays the role of the thermal vacuum 
energy. 
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The thermal excitations can then be described by a 
dispersion equation 



E 2 ,(k,T) = k 2 +m?(T), 



(14) 



where the subscript i = g,q labels gluons and quarks, 
respectively. The thermal masses rrn{T) are obtained 
from the self-energy of the corresponding particle, eval- 
uated at thermal momenta k ~ T. Eq. (14) is valid as 
long as this self-energy is only weakly momentum depen- 
dent in the relevant kinematic region. In addition, for a 
quasi-particle to be a meaningful concept, we require its 
thermal width to be small. The quasi-particle mass, for 
instance of gluons, is then 



m g {T)=T 



No 
6 



N 



+ ^l~g(T,N c ,N f ) (15) 



with the effective coupling specified as 

^• N "^ = 7m=w f ( il+ ^-W' (16) 

N c and Nf stand for the number of colors and flavors, re- 
spectively. The functional dependence of m g (T) on T is 
based on the conjecture that the phase transition is sec- 
ond order or weakly first order which suggests an almost 
power-like behavior m ~ (T — T c ) 7 with some pseudo- 
critical exponent 7 > 0. Setting g = 9.4, S = 1CP 6 
and 7 = 0.1, the effective mass, as given in eq. (15), ap- 
proaches the HTL result at high temperatures. The ther- 
mal quark mass is modeled similarly. Note that in con- 
trast to previous quasi-particle models extrapolated from 
HTL calculations [35] , the thermal masses used here drop 
as T c is approached from above. Obviously for T T c 
the perturbative limit of m g (T) and m q (T) is recovered. 
The density of the QGP takes the form 

i^C r d 3 k r d 3 k 

"( T ) = E 7(^)1 U(k,T) + Jj^s f B (k,T), (17) 

where the distribution functions for quarks and gluons 
are given by 

f qt , g (k,T) = i> q , g C(T) {exp [E qtt9 (k,T)/T] ± l} 

(18) 

being v q = 12 and v g — 16 the number of spin, flavor and 
color polarizations of quarks (Nf = 3) and gluons, while 



C(T) = Co 



(1 + S c ) 



(19) 



is the confinement factor. In order to reproduce lattice 
QCD thermodynamics for two light quarks and one heavy 
quark, the parameters take the values Co = 1.16, S c — 
0.02 and j c = 0.29. Pressure, energy density and entropy 
density of the QGP follow accordingly. 
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FIG. 2: Particle densities as a function of the temperature 
T in the hadronic and partonic phases. Here «sps denotes 
the initial quasi-particle gluon density at SPS and n^ HIC the 
corresponding value at RHIC. The vertical line indicates the 
value of the critical temperature T c = 170 MeV. 



The driving force of the transition at T c , the con- 
finement process, is physics that we must include phe- 
nomenologically. Below T c , the relevant degrees of free- 
dom are pions and heavier hadrons. Approaching T c from 
below, deconfinement sets in and quark and gluon quasi- 
particles are liberated. Conversely, when approaching the 
phase transition from above, the number of thermally 
active degrees of freedom is reduced due to the onset 
of confinement. As T comes closer to T c , an increas- 
ing number of quasi-particles gets trapped into hadrons 
(or glueballs). All confinement does on a large scale is 
to effectively reduce the number of thermally active de- 
grees of freedom as the temperature is lowered. This 
effect can be included in the quasi-particle picture by 
modifying the distribution functions by a temperature- 
dependent factor C(T) as given in eq. (19). Its explicit 
form is obtained by calculating the entropy density of the 
QGP with masses as in eq. (15). Dividing the lattice en- 
tropy density by this result yields C(T). Having fixed all 
parameters by extrapolating available lattice results to 
realistic zero temperature quark masses, the EoS for two 
light and a strange quark can be constructed for T > T c . 

Below the critical temperature, we employ the EoS of 
a non- interacting hadronic resonance gas, including all 
particles up to 1.6 GcV mass. A pion chemical poten- 
tial is later introduced as a mean to populate the pion 
phase-space occupation and reproduce the corresponding 
observed yield. Although this prescription is oversimpli- 
fied, any improvement has no bearing on the final result 
of this work. A matching with the high temperature 
part of the EoS is achieved by smoothly switching off the 
hadronic degrees of freedom. 

It is instructive to plot the total particle density given 
by eq. (17) as a function of temperature as shown in 
Fig. 2. The confinement factor C(T) accounts for the 
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release of the constituents of the hadrons as T grows. 
Since it does not jump abruptly to one as T c is reached, 
some hadronic clusters and hence non-zero hadron den- 
sities persist slightly above the critical temperature, as 
well as quasi-particles persist slightly below T c . We have 
indicated the values of particle densities at the beginning 
and the end of the time evolution of the fireball, as will be 
discussed in the next section. Comparing the numbers in- 
volved, it is obvious that the number density in the QGP 
is always at least an order of magnitude larger than the 
hadronic one (note the logarithmic scale). This has a 
crucial impact on *S? evolution in the produced medium. 

B. Time Evolution of the Fireball 

Having characterized the matter which constitutes a 
strongly interacting system in the thermodynamic limit, 
we now address the intriguing issue of how this system 
evolves in time after its production in a heavy-ion colli- 
sion. Our main assumption is the following: the fireball 
reaches local thermal (though not necessarily chemical) 
equilibrium within a time scale of the order of 1 fm/c. 
We then model the evolution dynamics by calculating the 
thermodynamic response of the hot and dense matter to 
the expansion of the total volume. Thermodynamical pa- 
rameters such as pressure p and energy density e in turn 
feed back in the expansion dynamics. 

We average all quantities over the fireball volume, leav- 
ing us with a spatially homogeneous system. The volume 
itself is taken to be cylindrically symmetric around the 
beam (z-)axis. In order to account for collective flow, 
we boost individual volume elements inside the fireball 
volume with velocities depending on their position. As 
flow velocities in longitudinal direction turn out to be 
close to the speed of light, we have to include the effects 
of time dilatation. On the other hand, we can neglect 
the additional time dilatation caused by transverse mo- 
tion, for which the velocity is typically v± <C v z . The 
thermodynamically relevant volume is then given by the 
collection of volume elements corresponding to the same 
proper time r. In order to characterize the volume ex- 
pansion within the given framework, we need first of all 
the expansion velocity in longitudinal direction as it ap- 
pears in the center of mass frame. Then we can compute 
the front position of the expanding cylinder as 

z(t) =v° z t+ I dt' f dt" a z (t") , (20) 

J to J to 

where a z (t) = c z p(t) / e(t) is the longitudinal acceleration 
and c z is a parameter. The time t starts running at to 
such that Zo = v z to is the initial longitudinal extension, 
v® being the initial longitudinal expansion velocity. The 
longitudinal position z(t) and t itself define a proper time 
curve t = yjt 2 — z 2 (t). Solving for i = t(r) one can 
construct z(t) — z(t). Then the position of the fireball 
front z(t) in the center of mass frame is translated into 



the longitudinal extension L(t) of the cylinder on the 
curve of constant proper time r. One obtains 




At the same proper time we define the transverse flow 
velocity and construct the transverse extension of the 
cylinder as a disc of radius 

R(t) = R + f dr' f dr" a ± {r") (22) 

Jto Jto 

where Ro is the initial overlap root mean square radius 
of the two colliding nuclei, a±(r) = cj_p(r)/e(r) is the 
transverse acceleration and c± is a parameter. With val- 
ues of L and R obtained at a given proper time, the 
3-dimensional volume is parametrized as 

V(t) = 2ttR 2 (t) L(t), (23) 

where the factor 2 arises from the relation Rbox = 
\/2R rms between the radius of a sharp cylinder and its 
root mean square value. 

The model is in principle fully constrained by fitting 
four parameters: the initial longitudinal velocity v®, the 
two constants c z , c± in a z , a± and the freeze-out proper 
time Tf. These are determined requiring 

R(Tf) = Rf, ™(Tf) = v f ±> ^(t(Tf))=v f z , T(r f ) = T f , 

(24) 

where Rf , v±, v( and Tf are extracted from an appropri- 
ate set of experimental data. In practice this is achieved 
only at SPS energy for central Pb + Pb collisions, where 
the freeze-out analysis [36] allows a complete determina- 
tion. We then assume entropy conservation during the 
expansion phase, fixing the entropy per baryon so from 
the number of produced particles per unit rapidity. Cal- 
culating the number of participant baryons as 

N p (b) = Jd 2 r T A (r)[l - cxp [- a ™T B (r)]} + (A <-> B), 

(25) 

where a™ = 30 mb at y/s = 17.3 GeV (SPS) and 
a™ = 42 mb at y/s = 200 GeV (RHIC), and scaling 
the entropy with the calculated number of participants, 
we find 5*0 = s A p , the total entropy. The entropy 
density at a given proper time is then determined by 
s(t) = So/V(t). Using the EoS given by the quasi- 
particle approach, thereby providing constraints from lat- 
tice QCD, we find T(s(t)) and also p(r) and e(r). The 
set of equations (20-23) is then solved iteratively, keeping 
the entropy constant. Finally we arrive at a thermody- 
namically self-consistent model for the fireball which is, 
by construction, able to reproduce the hadronic momen- 
tum spectra at freeze-out. This framework was shown 
to be a key ingredient for a successful description of the 
low- and intermediate-mass dilepton yields in SPS colli- 
sions (sec [18] where the model is also described in more 
detail) . 
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While the fireball evolution is fully constrained by data 
for central collisions at SPS, an extension towards differ- 
ent impact parameters and higher beam energies is not on 
the same firm ground. Unfortunately no detailed freeze- 
out analysis has been done so far for different centralities, 
even at SPS. One has to rely on sensible assumptions in 
order to construct the fireball evolution. 

Concerning the impact parameter dependence, from 
geometric considerations one can compute the overlap 
area and convert it, for simplicity, to a disc, while main- 
taining cylindrical symmetry but ignoring all effects of 
elliptic flow. This approximation should still be valid for 
not too large impact parameters. Note that the model 
is not applicable for very peripheral events, since the as- 
sumption of thermalization ceases to be valid. We then 
make the following two assumptions: first, the initial ve- 
locity w° of the fireball front is an important quantity 
which must change for more peripheral collisions. We 
consider the empirical fact that in pp collisions leading 
particles loose on average about one unit of rapidity. This 
is the limit we expect for very peripheral collisions. In 
order to account for this effect, we assume that the rapid- 
ity loss from incoming nuclei to the bulk of the produced 
matter scales with the number of binary collisions per 
participant (<~ 2.7 for central AB collisions, 1 for pp) and 
interpolate linearly the value of w° between these limits. 
Second, we assume that the freeze-out temperature Tf re- 
mains unchanged for different impact parameters. Since 
the total entropy scales with the number of participants 
N p , all parameters are fixed and the fireball evolution 
can be generalized to non-central collisions. 

The extension of the model to higher beam energies 
poses a more difficult task. The main reason lies in 
the transition from a baryon-rich to a baryon-poor sce- 
nario. While such an analysis is certainly possible with 
the amount of data collected by RHIC so far, it is a dif- 
ficult task and has not yet been carried out. For the 
time being, we aim at a qualitative description and re- 
frain from giving detailed predictions. Assuming that 
the entropy per participant scales with the total multi- 
plicity, we can calculate the initial entropy So. In order to 
account for the decreasing baryon density inside the fire- 
ball, we assume a rise in the freeze-out temperature from 
f 00 to f 30 MeV when going from SPS to RHIC energy. 
Similarly, we change the EoS in the hadronic part to in- 
corporate the decreasing pion chemical potential with in- 
creasing energy (for a more detailed discussion, see [18]). 
The thermalization time decreases from 1 fm/c at SPS 
to 0.6 fm/c at RHIC. From dN/dy particle distributions 
we extract the rapidity extent of hadron spectra to be of 
5.5 units, allowing v[ to be extracted. Moreover, exper- 
imental results indicate that transverse flow and the ge- 
ometrical freeze-out radius are virtually unchanged from 
SPS to RHIC energies [37], so we keep these two param- 
eters unchanged from the SPS case. Finally we adjust v® 
consistently with the assumed freeze-out temperature. 

The model parameters are eventually determined with 
the self-consistent procedure described above. Relevant 
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FIG. 3: Temperature T as function of proper time r for vari- 
ous collisions at SPS and RHIC energies. Curves correspond 
to impact parameters b = 0-1 fm (thicker curves), 2-3, 4-5, 
6-7, 8-9 and 10-11 fm (thinner curves). 



values of parameters for central collisions at SPS and 
RHIC are reported in Tab. I. One notices the significant 
QGP lifetime of about 8 fm, almost half of the total life- 
time of the fireball. Temperature T profiles as function 
of proper time r for collisions at different impact parame- 
ters at SPS and RHIC energies are plotted in Fig. 3. Note 
the delayed cooling occurring at the transition tempera- 
ture T c = 170 MeV. This is due to the softening of the 





SPS Pb+Pb 


RHIC Au+Au 




= 17.3 GeV 


y/s = 200 GeV 


TO 


1.0 [fm/c] 


0.6 [fm/c] 


T c 


8.0 


8.7 


T f 


18 


17 


To 


320 [MeV] 


390 [MeV] 


Tf 


100 


130 


Ro 


4.8 [fm] 


4.8 [fm] 


R c 


5.2 


6.0 


Rf 


9.0 


9.0 



TABLE I: Key properties of the fireball evolution deduced for 
central collisions at SPS and RHIC. The numbers given for 
the radii are root mean square values. 
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EoS, influencing the volume, and therefore temperature, 
expansion pattern via the acceleration a(r) oc p(r)/e(T) 
in eqs. (20) and (22). The presence of a longitudinal 
acceleration term also leads to interesting consequences. 
The smaller initial rapidity distribution of matter implies 
a larger initial energy density and hence a higher temper- 
ature of the fireball as compared to the Bjorken estimate 
for these quantities. 

Particle densities at the corresponding initial tempera- 
tures for central collisions at SPS and RHIC are indicated 
in Fig. 2 in the previous section. 



IV. CHARMONIUM EVOLUTION IN THE 
EXPANDING MATTER 




FIG. 4: Diagrams contributing to lowest order to ^ disso- 
ciation. The process involving a gluon quasi-particle (a) ex- 
hausts the amplitude since it already takes into account the 
one involving a quark quasi-particle (b). 



A. Interaction of \I> with the produced medium 

After having discussed the initial conditions for open 
and hidden charm production in AB collisions and the 
physics of the expanding medium, we now come to the 
central issue of the paper and examine the interaction be- 
tween charmonium and the produced medium. Because 
of the relatively small value of the J/ tp formation time 



/ 



{mj/^-m^') ^O.Sfm, 



(26) 



in the following we will assume that J/ ^>s are fully formed 
hadrons by the time the thcrmalizcd medium is pro- 
duced, and they will subsequently interact with its de- 
grees of freedom. Since the medium itself consists of a 
QGP for a significantly long time, we begin discussing 
how J/ip interacts with quark and gluon quasi-particles. 
As mentioned in the Introduction, we neglect the pos- 
sibility of static Debye screening, thereby focusing our 
analysis on the collisions of \P with the QGP degrees of 
freedom. We argue that a proper treatment of screening 
and collisions (together, in order to avoid double count- 
ing) is a complex dynamical problem involving an un- 
known time-scale for the modification of the cc binding 
potential. We therefore refrain from speculating on this 
issue, which deserves a separate study, and simply as- 
sume that "Js produced by the initial hard processes will 
propagate and collide with the QGP constituents. 

It is clear that collisions of with either quarks or 
gluons will lead to dissociation of the bound state. One 
expects that the two processes illustrated in Fig. (4) 
contribute to ^ dissociation in leading order. On the 
other hand, a quark can interact only via gluon exchange. 
Within the quasi-particle model, the process labelled (b) 
in the figure is effectively already included in the defini- 
tion of the temperature dependent gluon mass. Comput- 
ing both contributions would lead to an erroneous dou- 
ble counting. In other words, ^s only see gluonic quasi- 
particles in the plasma. As an exception, processes such 
as *c — ► gc and *c — > gc could be important. For exam- 
ple, the first one can proceed by dissociating — > cc and 
annihilating the c with an incoming c to form a gluon. 



Nevertheless such processes can be neglected when com- 
pared to gluon dissociation due to the much larger den- 
sity of the gluons themselves. 

We concentrate now on the process (a) in Fig. (4) for 
the reaction \f(p) + g(k) — ► c(qi) + c(q 2 ), and label mo- 
menta as indicated in parentheses. We then come to the 
problem of computing cross sections involving a relativis- 
tic bound state. In the present case one can argue that 
the cc system is, to first approximation, non-relativistic, 
greatly simplifying the treatment. Moreover, as was done 
originally by Bhanot and Peskin [19, 20] one can argue 
that the lowest-lying quarkonium levels can be approxi- 
mately described using the Coulomb part of the poten- 
tial. Then, with operator product expansion methods or 
more recent non-relativistic factorization techniques [38] , 
it is possible to obtain an analytic expression for the cross 
section. For the case of J/tj) the result is 



vj/tl>g(v) = Ao 



{uj/ej/^ - 1 ) 
(w/ej/^) 



3/2 



(27) 



where Aq = (2 n 7r/27) {m^tj/^) 1 ^ 2 - The cross section 
is a function of the gluon energy oj in the J / ip rest frame 
and involves the threshold energy ej/^, related to the 
binding energy ejy^, by the condition (p + k) 2 > Am 2 ., 
implying that e, 7/l/ , = £j rt + (e^) 2 /( 2m .w)- Tnc bind- 
ing energy is taken to be e j,^ — 800 MeV and the charm 
quark mass m c = 1.95 GeV as in [20], to fit the masses 
of first two charmonium levels (mju = 3.1 GeV and 
m^i =3.7 GeV). In an analogous fashion, it is possible 
to compute the cross section for tp' dissociation by gluons. 
Given the binding energy e°, = e^/4 (for a Coulombic 
system), and the cross section is 

n (,A ^KA (^/^' - 1 ) 3/2 (g/fg ~ 3) 2 . . 
a^ g {u) = 16A — . (28) 

Recently also the dissociation cross section for \c states 
was computed [39]. With a mass m Xc = 3.5 GeV and a 
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FIG. 5: Dissociation cross sections for the processes ^> g — > cc, 
where * = J/tp, ip' and Xc, as function of the gluon energy 
in the rest frame of >]/. 
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FIG. 6: Formation cross sections for the process cc — > tyg, 
where * = J/i/>, V>' and Xc, according to eq. (30), as function 
of the laboratory momentum in the rest frame of the target 
c. 



binding energy e Xr = 400 MeV one obtains 



(co/e 



1) 1/2 9(u,/e Xc ) 2 -20(u,/e Xc ) + 12 



(w/e x J 



(29) 

The cross sections calculated with eqs. (27), (28) and 
(29) are not necessarily reliable since the Coulomb po- 
tential approximation is used without justification. This 
problem is serious especially for the ip' and Xc states, 
which have much too large dissociation cross sections. 
Although the inclusion of higher states is in principle nec- 
essary in order to have a reliable description of the dis- 
sociation process, since ip' and Xc constitute about 40 % 
of the final J/ip yield, it is not possible at this stage to 
pursue further this approach. A way out could be to com- 
pute these cross section using a realistic potential for the 
cc system, but then questions arise for a rigorous separa- 
tion of perturbative and non-perturbative contributions 
and such a study falls beyond the scope of the present 
treatment. Within the adopted framework we choose to 
re-scale the tp' and Xc cross sections by a multiplicative 
factor k such that the obtained J ftp suppression pattern 
at SPS energies is in agreement with experiment. The 
required value is k = 0.2 and the resulting cross sections 
are plotted in Fig. (5). One notices the great difference 
in magnitude between the J /ip> case as compared to ip' 
and Xc even after re-scaling. We want to stress the fact 
that the value n = 0.2 is not necessarily large. In fact 
cross sections are proportional to the square of a matrix 
element. In the present case, since a gluon is a vector par- 
ticle, the matrix element consists in a derivative of the >F 
wave function in momentum space. A realistic confining 
potential would squeeze the spatial wave function more 
than the Coulomb potential does. In momentum space 
the wave function would then be more extended and flat- 
ter, resulting in a much lower cross section. Within the 
adopted framework we denote with o^ g ^ C c (s) the dis- 



sociation cross section, the dependent variable being the 
center of mass energy s, related to the gluon energy in 
the \t rest frame asw= k^p^/my = (s—p 2 — k 2 )/(2 m*). 

We now consider the possibility of cc coalescence in 
the QGP, a manner for production of 5* which has been 
recently considered by several authors. Processes such as 
gc — > and gc — > "J/c could in principle also contribute, 
due to the large number of gluons available. On the other 
hand they require the creation of a cc pair from the vac- 
uum and are therefore suppressed. In our formulation 
we then take into account only the fusion process. We 
do this by means of a cross section, applying detailed bal- 
ance to the reaction tyg <-> cc, and use the cross section 
calculated above for ^ dissociation by gluons. Choosing 
the zero momentum frame, flux factors are identical for 
the direct and reverse processes and simple kinematics 
guarantees that the relation 



4 (s-m%) 2 
3 s(s-4m 2 ) 



(30) 



holds. The factor 4/3 arises from counting the number 
of degrees of freedom (spin and color factors) in the two 
different channels. For the ^g system they are (1 x 3) x 
(8x2) = 48, while for the cc system (3 x 2) x (3 x 2) = 36. 
Since s = 2m c (E c +m c ) for the cc system, one can express 
eq. (30) as function of the c quark momentum pi ab in the 
rest frame of the c. The result for <7 ce ->* 9 is given in 
Fig. (6), where the formation cross sections for J/ip, tp' 
and Xc are shown. 

Coming now to the case of hadronic dissociation, we 
recall that many approaches have been developed in the 
literature to compute mesonic dissociation cross sections 
of ^. The first calculations were presented in [40]. Sub- 
sequently, more sophisticated treatments were developed 
(See [41] as an example). The resulting thermally aver- 
aged cross section amounts to a few mb for the process 
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"fp — > D.D, but computations considered p to be a par- 
ticle with infinite life-time, neglecting its broad spectral 
function. Accounting for this should considerably reduce 
the value of the calculated cross section. Significantly 
smaller values are obtained for ^tt — > DD. Despite the 
efforts applied to perform the mentioned calculations, a 
realistic time evolution of the fireball was never employed 
to obtain results comparable with experimental data. We 
stress that a proper and independent description of the 
produced medium is crucial if one wants to make use of 
charmonia as probes of dense matter. In our approach 
the medium evolution is constrained by the freeze-out 
analysis, leaving no room for adjustments. According 
to the calculation of the particle density as function of 
temperature, plotted in Fig. (2) in the previous section, 
we see that particle densities in the hadronic phase are 
almost two orders of magnitude lower than in the QGP 
phase. Note the logarithmic vertical scale. Unless one 
employs exceptionally large cross sections, it seems un- 
likely that hadronic dissociation can be at all relevant. 
A way out was proposed in [42], where it was argued 
that the QCD analogue of the Mott transition in an 
electron plasma may produce an onset behavior in ther- 
mally averaged cross section when the Mott temperature 
Tiviott ~ T c is reached. However, this framework has 
not yet been constrained by lattice EoS results. Given 
the inherent uncertainties and the above remarks on the 
smallness of particle densities in the hadronic phase, we 
will not consider the possibility of hadronic dissociation 
of*. 



Kinetic Description of Dissociation and 
Formation 



After having discussed and characterized the medium 
produced in a heavy-ion collision and specified how char- 
monium interacts by collisions with its constituents, we 
are ready to address the problem of how produced char- 
monia evolve and propagate through the medium until 
their final experimental observation. The natural frame- 
work in which to study the time evolution of * is that of 
kinetic theory. We use a semi-classical treatment, setting 
up a relativistic kinetic equation for the * phase-space 
distribution. The general problem of kinetic evolution 
of bound states in a strongly interacting medium was 
studied recently [43]. Here we make use of those results, 
adapting the treatment to our needs. We then assume 
that the kinetic equation describing the time evolution of 
the phase-space density of * is 

p^ fc(p, x) = C*(p, x) - C% (p, x) h(p, x) , (31) 

consisting of the drift term p^d^ on the l.h.s. and of 
a collision term on the r.h.s. The collision term is made 



of a dissociation (loss) part 

C%{p,x) = \ Jd<S> 3 (k,q u q 2 ) (2ir) 4 S\p+k- qi -q 2 ) 



X Wq,g^ C c fg(k,x) 



(32) 



describing * dissociation by quasi-particle gluons, stud- 
ied in [44, 45] to address the problem of charmonium 
suppression, and a formation (gain) part 

CfM = \ Jd$ 3 (k, qi ,q 2 ) (2n) i 5 i (p+k-q 1 -q 2 ) 

X Wcc^g fc( q i,x) f s ( q2 ,x) (33) 

describing W formation by cc fusion, introduced in a sim- 
plified manner in [22]. We stress here that eq. (31) stands 
actually for 3 equations, for J/ip, ip' an d Xc- 

In the equations above / c , f 5 and f g are the phase- 
space distributions of the degrees of freedom participat- 
ing in the collisions. Note that effects of Bose enhance- 
ment for gluons and Pauli blocking for c and c quarks are 
negligible. The space-time 4- vector is denoted by x, while 
kinematics is the same as in the previous section, that is 
p is the 4- momentum of *, k is that of the gluon, while 
qi and q2 arc those of the c and c quarks. The Lorentz 
invariant 3-body the phase-space integration measure is 



d$ 3 (k, qi,q 2 ) 



d 3 k 



d 3 



?1 



d 3 q 2 



(27r) 3 2^ fc (2tt) 3 2E 1 {2it) 3 2E 2 



(34) 



All particles are assumed to be on their mass shells. 

The transition probabilities W are averaged over the 
initial color and spin polarizations and summed over the 
final ones. Without the bar W indicates a transition 
probability already summed over both initial and final 
polarizations. Therefore W c s^^g = W c e^^g and 
W^g^ cS = W%g^ cE , with N cS = 36 and N 9g = 48. 
The transition probabilities satisfy detailed balance as 



W, 



(35) 



They are related to the cross sections discussed in the 
previous section, which can be expressed as 



°l (*) = -!- / d*2 (<Zi , 92 ) (2tt) 4 <5 4 (p+ k- q x - q 2 ) H* s _ 



for dissociation and 



1 



(36) 



4F C( 



a F( s ) = TEr- / d* 2 (p,*0 (27r) 4 S 4 (p+k- qi -q 2 ) W cS ^ g 

(37) 

for formation. Here the 2-body phase-space integration 
measure is 



d<$> 2 (p a ,Pb) = 



d 3 p a d 3 p b 
{2n) 3 2E a (2n) 3 2E b ' 



(38) 



while the flux factors are F cd = (p c ■ p d ) 2 — m 2 c vrp- d . In 
the non- relativistic approximation used in [38] , the tran- 
sition probabilities depend only on the center of mass 
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energy s. This allows to relate them to the cross sections 
as 

Vab^M = J^— W ^cd ■ (39) 

16tts F ab (s) 

It is then trivial to see that eqs. (35), (36) and (37) are 
consistent with eq. (30). 

Concerning the definitions of the phase-space distribu- 
tion of gluons, we adopt the expression given in eq. (18) 
of section III A. The time evolution of the gluon distri- 
bution is all contained in the proper time dependence of 
the temperature. For c quarks, as they are produced in 
the hard initial collision, we take their spectrum to be 
given by eq. (4). We assume that they do not interact 
significantly with the medium, approximately moving on 
straight lines according to the free streaming equation 



(40) 



The initial spatial distribution is assumed to be very nar- 
row. Introducing for convenience the space-time rapidity 
variable f] = 1/2 log[ (t + z)/(t — z) ], we define the phase- 
space density of charm quarks as 



fc{q,x) 



(2tt) 3 dN 



AB 



rm. 



dy c d 2 q A 



S(y c - V ) P± (r ± (T))e(T(T)-T c ) 

(41) 

and analogously for c quarks. The step function 9 is in- 
troduced to account for the fact that below T c c and c 
quarks have hadronized into D mesons and are no more 
available to coalesce into 'J. For simplicity in later nu- 
merical computations, the transverse position density is 
taken as a box of the same radius of the fireball as 



ttR 2 (t) 



9 (R(r) r ± (r)) , (42) 



where f_\_(r) = r j_ — (qj_/m € L ) t. The (^-function, arising 
from the assumption that c quarks are produced in a very 
narrow longitudinal region, strongly correlates y c and tj. 
This is a strong constraint on cc coalescence and we will 
discuss its implications later on. 

All the dynamics is now well defined since we have 
specified the distributions of the medium constituents 
and those of the bound state constituents, together with 
the transition probabilities Wy g ^ cS . The kinetic equa- 
tion can be solved exactly in closed form [43]. However, 
the medium evolution has been constructed by averag- 
ing over the whole interaction volume. It is therefore 
not possible to exploit the full phase-space information 
contained in eq. (31). We therefore simplify the kinetic 
equation as follows. 

We recall that the momentum spectrum at a given 
constant proper time can be computed from the phase- 
space distribution by means of the Cooper-Frye for- 
mula [46]. Integrating out also the transverse momenta, 
we obtain the rapidity distribution 

~d~y~^ = (2^f d2p± f d2r±dvTm * C0sh (y-V) h(p,x). 

(43) 



If we now re-express the operator p^d^ appearing in 
eq. (31) by means of the variables y and 77, we obtain 



p^d^ = m_L 



cosh(y — 77) 



d_ 



-sinh(y-r?)-^- 
t or] 



(44) 

It is now simple to see that differentiating with respect 
to r the * spectrum given with eq. (43) is equivalent 
to integrating the l.h.s. of eq. (31) in all phase-space 
variables except y. More precisely, we obtain 



d_dN* 
dr dy 



(r) 



(27T) 



7rP j d2pi ~ j d2rj - dl1 T P^ d IJ-fo(P' X )- 



(45) 

Integrating the r.h.s. of eq. (31) is a bit more compli- 
cated. We first discuss the dissociation term and make 
use of its definition in eq. (32) and of the definition of 
dissociation cross section given in eq. (36) to write 



/d 3 k 
( 27r )3g fc ^ °d ^ f g ( fc ' ,T ) ( p ' 



x . 



(46) 



If we now integrate C* over all phase-space variables 
except y we obtain 



/ 



d 2 p ± d 2 r ± dii t C%(p,x) fa{p 7 x) ~ 
d 3 k 



j^y, V*g ( s ) fg T ( T )) J d 2 ri_dri tE^ U (p, x) 

(47) 



dN q, 

- A D (r) _(r), 
where the dissociation rate is defined as 
A£to = Jj0pV 9g *Z(s) f g (k,T(r)) 



v = 

P± - <Pj 



(48) 



To obtain our result we have made certain approxima- 
tions. First of all we have suppressed the argument x in 
f g since it is redundant, the medium being uniform, and 
replaced it with the time dependent temperature. The 
first step to arrive at eq. (47) is obtained approximating 
the p± dependence of V^ g — g (s) / Ek and ao with 
its mean value. We neglect for simplicity the dependence 
of the mean value on centrality and choose (p±) = 2.0 
GeV, whose magnitude is in the range of experimental 
measurements, for all calculations. Our results are not 
very sensitive to the precise value of (p±), especially after 
folding them with the fireball time evolution. In the sec- 
ond step we assumed that is peaked around values of y 
close to 77, as for the c quark distribution in eq. (41), and 
set y = which is the relevant value for later comparing 
our final results with experiment. The dissociation rate, 
after converting it to a dissociation width by changing di- 
mensions from fm^ 1 to MeV, is plotted in Fig. 7. There 
the temperature dependence for J/tjj, tp' and Xc is shown. 
One notices that the J/tp width is smaller than the one 
of Xc for temperatures below ~ 300 MeV. J/ip becomes 
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FIG. 7: Dissociation widths of J/ip-, ip' and \c as function 
of the temperature. Arrows indicate initial temperatures for 
SPS and RHIC central collisions. 
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FIG. 8: Formation rates of J/ip, ip' and Xc computed at the 
initial evolution time as function of the collision energy. The 
values increase by more than three orders of magnitude from 
SPS to RHIC central collisions. 



broader at higher temperatures, while the Xc width stays 
virtually constant. Notice that the effect on ip' is negli- 
gible. The different behaviors are related to the different 
forms and thresholds of the dissociation cross sections, 
as illustrated in Fig. 5. 

Considering now the formation term on the r.h.s. of 
eq. (31), we integrate C* in the same way as for C* over 
the phase-space variables and define the formation rate 



A|(r) 



1 



(2tt) 3 



d 2 p± / d 2 r±dr] t Cf(p, x) 



(49) 



y = 



Its detailed evaluation proceeds without approximations 
but goes through several steps. We leave all details in 
the Appendix and here report the final result which is 

Af (t) = iy d 2 p ± d 2 q^d 2 q 2 Kaf(s)S c (y c ,qi)S c (y c ,q 2 ± ) 

where if is a dimensionless quantity given in eq. (A. 6), 
while S c (y, q±) is the c quark spectrum as given in eq. (4). 
The resulting formation rate A*p = A*(to) is plotted for 
J/ip, ip' and Xc in Fig. 8 at the initial evolution time, as 
function of the collision energy. The values increase by 
more than three orders of magnitude from SPS to RHIC 
central collisions, thereby suggesting that the mechanism 
of 'A 7 formation by cc coalescence might be important at 
high energies. 

With all the simplifications mentioned before, we have 
managed to reduce eq. (31) to a simple first order differ- 
ential equation for the rapidity distribution of $ as func- 
tion of proper time. In fact, putting together eqs. (45), 
(47) and (49) we obtain that 



d dN-q, -g, q, dNq, 

" v ) = x F (y> v - x d(v, t ) -j— ( t ) ' 



dr dy 



dy 



(51) 



whose solution, valid at y = 0, is obtained with few ele- 
mentary steps and provides the final \& rapidity distribu- 



tion at y = as 

dNj = (dN$ 
dy | dy 



exp 



Hdr' Xl(y,r') 

J Tn 



(52) 



+ 



f Tf dr' Xf(y,r') exp [- f\" \*(y 



Again, the solution found holds for the different char- 
monia J/ip, ip' and Xc- The initial condition dN^/dy 
is given by eq. (11) and the rates A* and A* are given 
respectively by eqs. (48) and (50). Note that they are 
implicitly dependent on impact parameter and collision 
energy. The structure of this solution is self-evident. The 
first term describes the dissociation of initially pro- 
duced in the hard collision, with the usual exponential 
suppression acting at all times from To to Tf, while the 
second term describes formation of ^s from cc in the 
QGP, from the initial time To up to an intermediate value 
t' and their subsequent suppression from t' to r/, inte- 
grated over all values of r'. This last term may become 
important as soon as the number of charmed quarks is 
large enough. This is expected to be the case in AB 
reactions as the collision energy increases, since the c 
quark spectrum is proportional to the mean number of 
nucleon-nucleon collision N co u > 1000 for large nuclei, 
and because the magnitude of the spectrum grows by 
nearly two orders of magnitude from, say, SPS to RHIC 
energy, as shown in Fig. 1. 

As a final remark we stress that leaving the y- 
integration open is important since we are interested in 
computing the value of the final \t rapidity distribution 
at mid-rapidity, and not the whole yield, for which ex- 
perimental results are available and with which we want 
to confront our model. 
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FIG. 9: Top: time evolution of dissociation rates A*(t) for 
J/tp, tp' and Xc, defined in eq. (48), computed for central 
collisions at SPS energy. Bottom: same as above for the 
integrated dissociation rates A*(r). 



C. Results 

Using the elements of the calculation as discussed in 
the previous sections, we can now compute the time de- 
pendence of dissociation and formation rates. The disso- 
ciation rate depends on the fireball temperature which, 
in turn, depends on time. How T depends in detail on r, 
was evaluated for different impact parameters and colli- 
sion energies in section III B. For central collisions at SPS 
energies (for RHIC conditions the result is qualitatively 
similar), the dissociation rates A*, defined in cq. (48) 
are plotted against r in Fig. 9 (top). For J/tp it drops 
fast, initially approximately as 1/r, for \c the rate drops 
more slowly, while for tp' it is maximum at hadronization. 
The integrated rates, denoted by A*, constitute the ex- 
ponents of the suppression factors in eq. (52). They are 
also plotted in Fig. 9 (bottom). All curves level off after 
hadronization in a similar way and the final value for Xc 
dominates the others. 

Analogously, we computed the formation rates A* , de- 
fined in eq. (50), which depend on time via the factor 1/r 
and a factor l/i? 2 (r), contained in the quantity K given 
with eq. (A. 6). It turns out that the coalescence mecha- 
nism is ineffective at SPS energies, therefore we plot the 
formation rate against time, for central RHIC collisions, 
in Fig. 10 (top). For all charmonia the formation rates 
drop faster with increasing r than the dissociation rates. 
They are negligible well before hadronization occurs. The 
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FIG. 10: Top: time evolution of formation rates A*(r) for 
J/tp, tp' and Xc, defined in eq. (50), computed for central 
collisions at RHIC energy. Bottom: same as above for the 
integrated formation rates A% (r) weighted by the suppression 
factor as in the r.h.s. of eq. (52). 



integrated formation rate is folded with the exponential 
suppression factor and computed as in the second term 
on the r.h.s. of eq. (52), denoting it by A*. The result 
is plotted also in Fig. 10 (bottom). Obviously all curves 
level off before hadronization, showing a dominance of cc 
coalescence into J/tp rather than in tp' or \c- Noteworthy 
is the tendency of the A*s to decrease once again at later 
times, when the A*s are by then small and dissociation 
tends to take over. This behavior is more pronounced for 
Xc, for which dissociation is still strong in the vicinity of 
hadronization. 

Dissociation and formation rates are then combined to- 
gether according to eq. (52) to give the observed ^ rapid- 
ity distribution. As the fireball evolution is constructed 
for different impact parameters and collision energies, 
we compute the corresponding rapidity distribution at 
y = 0. 

Since we are interested for the final J/tp yield, we need 
to take into account decays into the charmonium ground 
state of tp' and Xc- Experimentally is has been observed 
that measured J / tps come in a fraction of 60 % from di- 
rect production, while a fraction of 10 % comes from tp' 
decays and a fraction of 30 % is from Xc- Denoting these 
fractions as Wj/^p — 0.6, w^> = 0.1 and w Xa — 0.3, we 
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can write the final J/ip rapidity distribution as 



dN 



fin 



dy 



E 



w _ _t 
vvy — — 

dy 



(53) 



where the sum extends over * = J/ip,tp' ,\ c - Approxi- 
mating the weighted exponential for suppression of the 
initially produced J/ips, implicit in cq. (53) and appear- 
ing in eq. (52), as 



£ w* dN * 



dy 



exp 



dy 



E 



f'dr' Xt(y,r') 

J To 

Hdr' \%(y,T') 

J TO 



(54) 



W<a exp 



which is correct to leading order in the expansion of the 
exponential, we obtain the final form of the J/ip rapidity 
distribution which can be compared with experimental 
data. The required expression is 



dN 



fin 



(55) 



y)W*exp - Tdr' Xl(y,r') 

d y { d y v L ^ 

+ E W * Heir' \%(y,r') exp [- f T dr" \%(y,r") 

,T, -'"TO I Jt' 



We now look at the case of Pb+Pb collision at ^/s = 
17.4 GeV and compute the J/ip spectrum as function of 
the impact parameter b. Then we construct the J/ip/DY 
ratio 



Rj/->p/DY(b) = B^+ 



dNj^(b)/dy 



B„ 



dN DY (b)/dy 
dN J/4 ,(b)/dy 

— iV 



(56) 



dNg Y /dy N coll (b) 



dN J/4 ,(b)/dy 
dNf/dyN cM (by 



where we assume that the Drell-Yan spectrum in Pb + 
Pb collisions scales with the number of collisions N co u. 
The coefficient B^^- is the branching of J/ip into 
pairs. The overall normalization is fixed at No = 53.5 to 
match, at large values of 6, the experimental value in pp 
collisions. The J/ip/DY ratio is then plotted in Fig. 11 
(solid line) as function of the mean transverse energy 



E T (b) = e T N p (b). 



(57) 



The quantity et = 0.274 is the amount of produced trans- 
verse energy per participant. As reference, we also plot 
the curve obtained by considering only nuclear effects 
(dotted line) and neglecting the contribution of the pro- 
duced medium. The agreement with data is quite re- 
markable, in particular the slope of the curve, considering 
that fireball parameters have not at all been tuned to this 
particular observable. The result as such deserves some 
detailed comments. First of all, for simplicity we do not 
perform the usual convolution with the Et — b correlation 
function [13], therefore it is obvious that the curves end 
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FIG. 11: Result at SPS energy for the J/ip/DY ratio as func- 
tion of the transverse energy. The dotted curve, labelled 
U NUC" , includes only nuclear effects, while the solid line, 
labelled "QGP" , is the complete result including dissociation 
by collisions with gluonic quasi-particles. 



at Et — 110 GeV, which corresponds to the mean trans- 
verse energy at b = 0. To go beyond this point and be in 
agreement with the data it is necessary to include effects 
of fluctuations, which can be treated in a straightforward 
manner [47]. Second, we find that the formation mecha- 
nism of J / ip by cc coalescence is totally negligible at this 
energy, while suppression is caused exclusively by colli- 
sions with gluonic quasi-particles. Third, hadronic disso- 
ciation of J/ip seems to be ruled out, not because cross 
sections are small, rather because the hadronic number 
density is more than an order of magnitude lower than 
the partonic one, and only a fraction of a fm -3 already at 
hadronization. All this was clearly shown in Fig. 2. This 
last statement relies on the fact that the time evolution 
of the fireball is not constructed ad hoc in order to obtain 
the desired result, but a priori and with independent ex- 
perimental information (hadronic and low-mass dilepton 
spectra). Last, the low Et region is over-estimated in 
our result. This is mainly due to the too simple scaling 
with the number of participants, neglecting the fact that 
for central collisions the amount of produced transverse 
energy receives a contribution from the mean number of 
collisions. In any case, the fireball description is not ex- 
pected to be valid at large impact parameters. 

As we now increase the collision energy, staying at im- 
pact parameter b = 0, we observe an interesting feature. 
Since the amount of charmed quarks grows substantially 
with increasing energy, as was indicated in Fig. 1 , we ex- 
pect that the formation of J/ips via cc coalescence can 
become important and at least comparable to the pri- 
mordial one. To study this we construct the suppression 
function 



q M dNj^{s,b = 0)/dy 
*J/^ S ) dN^(s)/dy N coll (b = 0) 



(58) 
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FIG. 12: Suppression factor Sj/^ as function of the collision 
energy. The thin solid line labelled "NUC" is computed with 
nuclear effects alone. The thin dotted line labelled "cjf = 0" 
is the result with only J/tp dissociation. The thick solid line 
labelled "Ay = 0" is the full calculation. The thick dotted and 
dot-dashed lines, labelled "Ay — 2" and "Ay — 4" respectively, 
are obtained by relaxing the assumption of complete y c — ij 
correlation, as explained in the text. 



and look at its dependence on s. Note that in absence of 



nuclear or medium effects one should have S 



J/V 



1. In- 



dications of the onset of the formation mechanism clearly 
appear from our calculations, which we plot in Fig. 12. 
The thin solid line labelled "NUC" represents the sup- 
pression function when only nuclear effects are included. 
According to eq. (13) the effective absorption cross sec- 
tion is a$ uc {yfs = 200 GeV) = 8.2 mb. The thin dotted 
line labelled u a F — 0" is the result with only * disso- 
ciation, without the formation mechanism. The thick 
solid line labelled "Ay = 0" is the full calculation. One 
clearly observes, especially at the full RHIC energy, a 
substantial contribution of about 30 % of the total from 
the formation mechanism. The thick dotted and dashed 
lines, labelled "Ay = 2" and "Ay = 4" are obtained by 
relaxing the assumption of complete y c — r\ correlation 
realized by the 5- function in eq. (41). This is done ap- 
proximately by smearing the latter with a Gaussian of 
width Ay. The values chosen are Ay = 0, the baseline 
value, and Ay = 2,4. More details are given at the end 
of the Appendix. Smearing the y c — ry correlation allows 
the formation mechanism to become more efficient, lead- 
ing to a J/tp yield of comparable magnitude with the 
primordial one. The net effect is a slowly rising function 
of y/s, but in any case we still have suppression, even at 
RHIC energy. 

We now proceed to examine central collisions at RHIC 
energy. In order to compare with preliminary data from 
the PHENIX experiment [31], obtained at mid-rapidity 
from ^ decays into e + e~ pairs, we construct the quantity 
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FIG. 13: Result at RHIC energy for the J/tp yield, scaled 
by the number of binary collisions, as function of the number 
of participants. Different line types for the curves and their 
labels are described in Fig. 12. Filled squares are preliminary 
PHENIX data. 



and plot it in Fig. 13 as function of the number of par- 
ticipants. The quantity B e + e - is the branching of J/tp 
into e + e~ pairs, while the overall normalization is fixed 
at N^(N t 



N*^(b) = B e+e 



dN J/4 ,(b)/dy 
N coll (b) 



(59) 



0) = 0.15. We label all curves consis- 
tently with Fig. 12. Although a comparison with the 
data is, at present, premature, we see that our results 
lie within experimental errors. Again, the contribution 
from cc coalescence into J/tp is significant, although not 
dramatic. All the obtained curves are monotonically de- 
creasing and we do not find any inversion of this tendency 
at any large value of N p . In other words, we do not find a 
net enhancement of J/tp, but this result needs to be con- 
firmed by more accurate calculations. The reasons are 
twofold: first, although we have a rough picture of what 
medium to expect at RHIC, the corresponding fireball is 
presently not yet determined at the same level of quality 
as with existing SPS data. Second, shadowing effects for 
J/tp production have been taken into account only very 
crudely. For open charm production also recent calcu- 
lations [28] show a non negligible amount of shadowing. 
Indeed, the expectation at RHIC energy is a reduction of 
dN% B /dy by about 30% at mid-rapidity. The resulting 
reduction factor of 0.7 enters squared in the formation 
rate given on the right of eq. (50), therefore reducing by 
about half the enhancement effect. 



D. Remarks on the QGP at SPS 

Consider once again the solution of the kinetic equation 
given as eq. (52) and neglect the formation term. Taking 
the time evolution of the medium density in the simple 
form ti(t) = uqTq/t and assuming that the average dis- 
sociation cross section is time independent (denoted by 
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(7d), wc obtain the simple analytical solution 

o\d n r log(n / nj ) , (60) 



dy 



dy 



exp 



where nj = n{rf). We examine the exponential suppres- 
sion factor for parameter values corresponding to the sit- 
uation in central collisions at the SPS. We assume that 
the produced QGP is thermalized so that the gluon den- 
sity is given, in the ideal gas limit, by 



n(T) 



ffC(3) 



(61) 



The final density is fixed by eq. (61) at the critical tem- 
perature T c = 170 MeV, yielding n f = 1.25 fur 3 . Tak- 
ing <jd = 1 mb and tq = 1 fm, consistently with what 
has been computed earlier in the previous sections, one 
needs hq = 5 fm~ 3 in order to have a suppression factor 
of about 0.5, which is precisely the amount of anoma- 
lous suppression required by the NA50 data at large Et- 
Note that the initial density required in this scenario cor- 
responds to T = 270 MeV. 

We now compare the above estimate with the results 
given by a popular model [48], describing J/tp suppres- 
sion by means of collisions with comovers, usually con- 
sidered to belong to the class of "hadronic models" . Its 
essential element is the suppression function 

S(b, r) = exp [-a co N c y °{b, f) log( N™(b, r)/N f )] , (62) 

where a co = 1 mb represents the interaction with co- 
movers, Ny°(b, s) is the 2-dimcnsional rapidity density of 
hadrons as function of the transverse coordinate r for a 
given impact parameter b, computed in the Dual Parton 
Model [49], and Nf = 1.15 fm~ 2 . Averaging S by means 
of the "production probability" TaTb and neglecting nu- 
clear effects, we obtain 



:s)(b) 



exp 



a co N y (b)log(N y {b)/Nf) , (63) 



where the average rapidity density of comovers is com- 
puted as 



Ny(b) -- 

If we now identify 



d 2 rT A (r)T B (r) N c y °{b,r) 



I- 



(64) 



d z r T A (r)T B (f) 



&d = <? co , n t = N y (b = 0) and n f t = Nf , (65) 

we immediately recover the suppression factor in eq. (60). 
Computing numerically eq. (64) at b = 0, we obtain 
iVj, = 6 fm~ 2 , which should be compared with the prod- 
uct no To = 5 fm~ 2 . The numerical agreement is quite 
striking. To show that this is not a mere coincidence, 
we can compare the final values for the densities. Nu- 
merically this means to confront Nf = 1.15 fm -2 with 



nfTo — 1.25 fm~ 2 determined above: again, a remark- 
able agreement. It is clear that the initial density of 
"hadrons" as computed with the Dual Parton Model is 
far too high for the medium to be "hadronic", and the 
particle density Nf is not to be understood as the freeze- 
out value but rather to be taken at hadronization. It is 
therefore suggestive to interpret the comover model as 
actually describing anomalous J/tp suppression by the 
QGP. 



V. CONCLUSIONS 

We have calculated J/tp production from an expanding 
fireball created in rclativistic heavy-ion collisions over a 
wide range of centralities and beam energies, from SPS at 
y/s = 17.3 GeV (Ei ab = 158 GeV) to RHIC at y/7s = 200 
GeV. In this summary we review our basic assumptions 
put them into perspective. 

The produced medium was described assuming QGP 
formation and thcrmalization. A phcnomcnological 
quasi-particle model for quarks and gluons, in accor- 
dance with lattice QCD thermodynamics, was applied to 
model the partonic phase. This provided a realistic EoS 
which we then used to drive the expansion dynamics of 
the medium by means of a fireball model, characterized 
by time dependent temperature and volume. This con- 
struction had already been proven successful in describ- 
ing low mass dilepton spectra at SPS. We want to stress 
again that the setup of this fireball (initial and freeze- 
out temperature, density etc.) had been fixed by inde- 
pendent observables and ensures consistency with a mul- 
titude of hadronic measurements. In contrast to previ- 
ous approaches, we have therefore largely eliminated the 
medium evolution as an adjustable parameter in the J/tp 
description. Nevertheless, this framework is still simpli- 
fied and is probably insufficient to account for phase- 
space correlations of the medium, which have been av- 
eraged out. In a future attempt to study rapidity and 
p^-dependent properties of J/tp, this is a problem to be 
taken into account. 

We then set up a kinetic equation with a collision term 
incorporating both gain and loss terms. These stand for 
\]/ formation due to coalescence of cc quark pairs and H! 
dissociation due to collisions with gluon quasi-particles, 
respectively. The elementary process ^g — > cc was mod- 
eled by a simplified dissociation cross section, approxi- 
mating $ as a Coulomb bound state. The corresponding 
back reaction cc — > ^g could then be obtained by de- 
tailed balance. Here a few important remarks are neces- 
sary. First of all, a crucial ingredient for the construction 
of the solution of the kinetic equation is a reliable initial 
condition. In the present work we have simplified this 
aspect in order to focus on the QGP effects, but it is 
possible to improve on this point, especially in what con- 
cerns the description of nuclear effects, both for open and 
for hidden charm. New data collected at RHIC in recent 
deuteron+Au collisions will certainly help constraining 
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the model dependence of such calculations. Second, the 
description of medium effects on \1/ was developed ne- 
glecting screening of the cc potential. This is an assump- 
tion which should be checked by a comprehensive study 
of the ^ self energy in a QGP, including, at the same 
time, screening and collision processes. Third, the use of 
the Coulomb part of the potential to describe cc bound 
states is only a rough approximation. In fact, cross sec- 
tions involving ip' and \c had to be rescaled in order 
to achieve agreement with the data. A proper descrip- 
tion of bound states with a realistic potential, including 
a confining and possibly screened part, is then necessary. 
Despite such simplifications, the description as a whole 
seems to be coherent and provides a framework for future 
developments. 

With proper averaging of the kinetic equation over the 
spatial extent of the fireball, we then found a simple solu- 
tion, which allowed direct comparison with experiment. 
At SPS energy we were able to describe the suppression 
effect in the data, without the need to invoke hadronic 
comovers. These results support the hypothesis that the 
QGP is actually produced, at a transient stage, in Pbj-Pb 
collisions at y/s — 17.3 GeV. It is important to note that, 
within our outlined approach, a purely hadronic frame- 



work would not be successful in describing existing data. 
Moreover, since the hadronic phase exists only at moder- 
ately low particle densities, it has no bearing on iff evo- 
lution. We also considered extrapolations up to RHIC 
energies where, despite the more extreme conditions as 
compared to SPS, a sizable fraction of primordial J ftps 
still survives. Although a clear trend towards more copi- 
ous J/ip production in by cc coalescence was found, no 
net J ftp enhancement was present in the end. This result 
needs to be confirmed by a dedicated study of space-time 
and momentum r\ — y c rapidity correlations for charm 
quarks, since they have a non-negligible effect on the fi- 
nal results. 
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APPENDIX: EVALUATION OF THE COALESCENCE RATE 



We evaluate here the coalescence rate. From the definition of the formation term in eq. (33), substituting in it 
eqs. (34), (39) and (41) we obtain 



C F (p, x) = ^5- / ^ ~iF~ 5 *(P + fc - 9i - 92) S c (yi,q{) S c (y 2 , q 2 ± ) 5{yi - rf) S(y 2 - tj) 



6(R(r) \f ± - (Si/mi) r|) 6(R(t) \f ± tf ± /m 2 ± )r\) . 



' J E k E 1 E 2 



(A.l) 



Integrating over all phase space variables, except for y, we define the formation rate as 
= (2^)3 Jd 2 p±Jd 2 r ± di 1 TC F {p 7 x) 

= 4^: J rf2fc -L A± d 2 q 2 ± J dy k dm dy 2 5{Ey + E k - E x - E 2 ) 5(p z + k z - q\ - q 2 z ) 5{ Vl - y 2 ) 

S (s - 4ml) 1 \ n I 1 \ o / 2 \ 

X m \m\{ s - m %) ^ (S) S c (y^)S c (y 2 ,q ± ) 

X f d2r± r I , l2 fl(^)-l^-|/ 2 i) 0(R(T)-\r ± + Z/2\), (A.2) 
J VkR (t)J 

where k ± = q\ + q\ — p± and £ = {q\/m\_ — (f^/m 2 ^) r. The integral over rapidities can be manipulated to give 



j dy k dyi dy 2 (5(m* + k cosh y k — m]_ cosh y\ — m\ cosh y 2 ) 
X (5(fcsinh y k — m]_ sinhj/i — m\ sinhy 2 ) $(yi — Mi) 
J d Vl dy 2 ^(yi - y c ) S(y 2 - y c ) 6{fi 2 N - fi 2 D ) 



V(^) 2 - U%T 2 



(A.3) 



18 



where p 2 N = (m*) 2 + (m]_ + ml) 2 - k\, n% = 2mf {m\ + m 2 ) and y c = log (p 2 N /p 2 D + ^/(Mat/m?)) 2 - ■ Thc 

integral over the transverse position variable can be computed analytically and it turns out to depend only on the 
variable x = £/2R(t). The result of the integration is 



Q(x) 



2 j 

1 \x\J 1 — x 2 + arcsinx 



itR 2 (t) 

In this way the formation rate can be brought to thc form 

,2 



(l-x). (A.4) 



1 f i ( <t — Am 2 ) 

X F (r) = - d 2 k ± d 2 qld 2 ql \ (S c> 2 , a F (s) S c (y c , q\_) S c (y c , q 2 ± ) (A.5) 

t J it mj_mj_ [s — m^) 

Hp 2 n - A)- 



Q( X ) a i . 2 



Defining the quantity 



VK) 2 - 0%) 



K _ s (s~ 4m 2 c ) Q(x) 2 2 



we finally arrive at the expression given in eq. (50). 

If we now relax the assumption that the c quark distribution contains the factor S(y c — rj), replacing it with a 
Gaussian g(y c — rj) = (nAy 2 )^ 1 ^ 2 cxp[— (y c — rj) 2 /Ay 2 }, the result of the 77-integral, which is one of those to be done 
to obtain eq. (A. 2), is 

J d V g( yi - rf) g(y 2 - v) = (lit Ay 2 )- 1 ' 2 cxp [-(y c - V ) 2 /2Ay 2 ] . (A.7) 

Defining p = y\ — yi and y c — (y\ + y-i)j2 and expanding the l.h.s. of eq. (A. 3) to order p 2 7 we approximate the 
various terms typically as 

m]_ cosh ?/i + m\ cosh 7/2 — (m]_ + m\) cosh y c (1 + p 2 /8) and (A. 8) 

m]_ sinh y\ + m 2 sinh 7/2 — (m±_ + rn\) sinhj/ c (1 + p 2 /8). (A. 9) 

We have assumed that m\ — m\ ~ for the values of q]_ and q\ which dominate the final phase-space integration. 
Substituting p 2 with the mean value Ay 2 , the final result is reduced to the substitution 

mi+mi— >(mi+mi) (1 + Ay 2 /8) (A.10) 

in p, 2 N . This amounts to open up phase-space for the transverse integrations since there are more values of transverse 
momenta which satisfy the constraint 9(p 2 N — p 2 D ). 
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